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Abstract 

Coexistence of individuals with different species or phenotypes is often found in nature in spite 
of competition between them. Stable coexistence of multiple types of individuals have implications 
for maintenance of ecological biodiversity and emergence of altruism in society, to name a few. 
Various mechanisms of coexistence including spatial structure of populations, heterogeneous indi- 
viduals, and heterogeneous environments, have been proposed. In reality, individuals disperse and 
interact on complex networks. We examine how heterogeneous degree distributions of networks 
influence coexistence, focusing on models of cyclically competing species. We show analytically 
and numerically that heterogeneity in degree distributions promotes stable coexistence. 



PACS numbers: 89.75.Fb, 87.23.Cc, 89.75.Hc 



I. INTRODUCTION 



How to maintain or prevent coexistence of competing multiple types of individuals is a 
key issue in various areas. For example, coexistence of multiple species in ecological habitats 
taplies stab.e biodiversity reaped in nature HQ. Coexistence of multiple types of p.ayer S 
in evolutionary games implies survival of altruistic players in the sea of selfish players [3J. 
Coexistence of disease-free and infected individuals implies an endemic state that should be 
suppressed usually [l]. 

Mechanisms of coexistence have been a central theoretical question because complexity of 



:ity or 



a population state (i.e. coexistence) and stability are often contradicting requirements 
Q]. Coexistence in population dynamics has been explained by, for example, nonequilibrium- 



state interpretation, habitat subdivision, heterogeneity in species such as heterogeneous 
dispersal speeds, and heterogeneity in environments p], 0, Q, fj- Spatial structure such 
as the square lattice also limits diffusion and promotes coexistence. In this case, each 
species is clustered in different regions of the lattice 0, [h], H, 12] . However, real- world 



interaction quite often occurs on contact networks of individuals that are more complex 
than the square lattice. Most real networks have the small-world and scale-free properties 
(e.g. Refs. Q])- The small- world property is equivalent to the combination of small 
average distance between vertices and large clustering, or abundance of densely connected 
small subgraphs such as triangles. A scale-free network is defined by a degree distribution, 
or the distribution of the number of contacts (edges) that each vertex has, which follows 
a power law, pk oc A;" 7 . Here pk is the probability that a vertex has degree k. The scale- 
free property may be too idealistic to describe contact networks underlying real population 
dynamics. Even so, it seems likely that different patches or individuals are endowed with 
different connectivity to others. 

In terms of networks, some known mechanisms of coexistence benefit from the regular 
lattices, the one-dimensional continuous line, the two-dimensional continuous plane, and the 
complete graph (mean field situation), in which all the vertices are considered to share the 
same degree (see the papers cited above and references therein). In random graphs, which 
are sometimes used in this context the vertex degree obeys the Poisson distribution. 
However, the real degree distribution may be even broader. 

Here we investigate how possibility of coexistence is affected by heterogeneous degree 
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distributions of contact networks, not by heterogeneous environments other than network- 
based ones or heterogeneous individuals. Among various competitive relationships among 
different phenotypes, we focus on cyclic competition of three species, which is a minimal 
case. 

Cyclic competition is actually abundant in nature. For example, tropical marine ecosys- 
tems 15] and vertebrate communities in high-arctic areas 16] include cyclic dominance 

n 

relations composed of a couple of organisms (also see Ref. [7J). Real microbial communities 



of Escherichia coli 



171 ] and color polymorphisms of natural lizards 18| also have cyclically 



dominating three phenotypes and show alternating wax-and-wane population densities. In 
evolutionary games, the public-good game with volunteering, namely, the choice of not 
joining the game, results in cyclic competition [19(. The susceptible-infected-recovered- 
susceptible model of epidemiology and models with additional types of states also include 
cyclic competition [4, 2(|. We focus on two specific predator-prey models of such cyclic 
interaction, that is, the standard rock-scissors-paper (RSP) model 7J and the May-Leonard 
(ML) model 2l|]. These models have neutrally stable or unstable coexistence solutions in 
well-mixed populations. Therefore, in a finite population, population dynamics are eventu- 
al. 



ally trapped by an absorbing state corresponding to the dominance of one species [?], 
We show that heterogeneous degree distributions stabilize coexistence of multiple types of 
individuals placed on networks. 



II. ROCK-SCISSORS-PAPER DYNAMICS ON NETWORKS 



A. Model 



As a minimal model of cyclic competition, we consider the standard RSP dynamics on 
networks with heterogeneous contact rates. There are three species, which we call states, 
represented by rock, scissors, and paper; rock beats scissors, scissors beat paper, and paper 
beats rock. Each vertex takes state 0, 1, or 2. A pair of vertices may be connected by an 
edge. The degree A; of a vertex is the number of edges, or the number of contacts with other 
vertices. State 1 outcompetes state by invading onto each neighboring state with state 
at a rate of A. In other words, a vertex with state changes its state to 1 at a rate of Ani, 
where n, is the number of vertices with state i in the neighborhood. Similarly, 1 (2) turns 



3 



into 2 (0) at a rate of pn 2 (n ). In an ecological context, we are considering the limit that 
dispersion rates (= 1, A, p) are much larger than the natural death rate. We consider the 
influence of death rates later with the ML model. 

For a perfectly mixed population, the mean field theory tells that there are an ensemble 
of neutrally stable limit cycles surrounding a neutrally stable equilibrium corresponding to 
coexistence of the three states. Therefore, the coexistence solution is practically unstable 
in finite populations. The RSP dynamics with spatial structure, such as the square lattice, 
accommodate many states each of which is clustered in different loci jlO . 
are interested in a network mechanism that may enable stable coexistence. 
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12] . Here we 



B. Equilibrium 

With dispersed degrees, vertices with different degrees obey different state-transition 
dynamics. Let us denote by pi jk the probability that a vertex with degree k takes state i 
(= 0,1,2). The probability that a vertex adjacent to an arbitrary vertex takes state % is 
denoted by 0j. This probability does not generally agree with p^ k or its average over all the 
vertices. This is because a vertex with more edges is more likely to be selected as a neighbor. 
In fact, a neighbor has degree k with probability kpk/ (k), where pk is the probability that 
a vertex has degree k and (k) = kpk is the mean degree giving normalization. Therefore, 
®* = J2k^PkPi,k/ (k) [23]. Because each vertex is occupied by one of the three species, 
namely, p ,fe — 1 — Pi,k — P2,k and O = 1 — Oi — 02, it suffices to consider the density of 
state 1 and that of state 2. Noting that the expected number of state-i neighbors of a vertex 
with degree k is equal to kQi, we derive 

p hk = A(l - p ljk - p 2 ,k)kQ 1 - Lipi, k kQ 2 , (1) 
p2,k = PPi,kkQ 2 - P2,kk(l - 0i - 2 ). (2) 

For example, the first term in Eq. (CQ) corresponds to the invasion of state 1 onto vertices 
with state 0. In the equilibrium, we have 

I Plk \ = A0| / 1 - e; - 65 

V Plk ) ( Ae i + - i - 2) + ^ ^ 
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The coexistence solution to Eq. (J3J) is given by 

i / 1 \ 

(4) 





20] 



'2 

for any k. The degree distribution does not affect the equilibrium population densities 
C. Stability of coexistence equilibrium 

When pk — $k,(k}j each vertex has the same degree equal to the mean (k). This case 
corresponds to well-mixed populations. Then the coexistence is neutrally stable (e.g. Refs. 
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21] ), which under 



13, 



ies experimental and natural ecosystems showing oscillatory popula- 



18] ■ Equation ([2]) indicates that the oscillation period is proportional 



tion dynamics 
to l/(k). 

However, the stability of the coexistence and realized dynamics are considerably influ- 
enced by networks. To see this, let us consider a two-point degree distribution given by 
Pk = p8k,ki + (1 ~ p)dk,k 2 - On average, a total of np vertices have degree k\ and n(l — p) 
vertices have degree k 2 . Equations and §2§ for a network with the two-point degree 
distribution define a four- dimensional dynamical system. We set A = /i = 1 for simplicity, 
although generalization to other A and fi is straightforward. The characteristic equation 
evaluated at the coexistence equilibrium [Eq. (TJ])] is represented by 

4 3hk 2 , fk 1 k 2 (k 1 + k 2 - (k)) (k 2 ) 2 \ 2 (k 2 ) , , l2l2 

x + H x +?j { V + w) + w i2X + 12 = ' (5) 

where (k) = kpk = pkx + ll—p)^ and (k 2 ) = k 2 Pk — pk 2 + (l—p)k 2 . When k\ — k 2 — k, 
we turn back to the ordinary mean field case with neutrally stable oscillations: x = \/3ki, 
(— 3 ± \^3)k/2. More generally, the Routh-Hurwitz criteria for Eq. (0) is 

\Hi\ = 1, 

\H 2 \ = pk 2 ((k) - h/2) 2 + (1 - p)k\ ((k) - k 2 /2f + Zklk 2 2 /A > 0, 

|# 3 | = 81k 2 1 k 2 2 p(l-p)(k 2 -k l ) 2 ^(k 2 ) 2 + 2k l k 2 (k) 2 )/(k)\ 

\H 4 \ = 9k 2 kl\H 3 \, (6) 

where \Hi\ is the zth principal minor. The coexistence solution is stable when \H 3 \, \H 4 \ > 0, 
that is, k\ k 2 and p^O, 1. Dispersed contact rates stabilize coexistence. 



D. Numerical results 



We resort to numerical simulations to examine more general networks and to be more 
quantitative about the effects of degree dispersion. We compare different types of networks 
with n = 5000 vertices and (k) = 10. The regular (R) random graph corresponding to the 
ordinary mean field case is generated by the configuration model 14| with pk = S^/k)- This is 
a type of random graph in which every vertex has the same degree (k). We also use the Erdos- 
Renyi (ER) random graph, which has the Poisson degree distribution pk = (k) k /k\, 
and the Barabasi- Albert (BA) scale-free network with the parameter m = (k) /2 = 5, which 
yields pk oc A; -3 (k > m) and pu = (k < m) [131 ] . 

Typical population dynamics for these networks are shown in Figs. [H(a)-[I](c). On the 
R random graph, the coexistence solution is neutrally stable in theory. Combined with a 
finite-size effect, the amplitude of the dynamical population density becomes progressively 
large in an oscillatory fashion. Eventually, one state dominates the whole network in an 
early stage [Fig. [1(a)]. On the ER random graph [Fig. [1(b)] and the BA model [Fig. DD(c)], 
coexistence occurs owing to the distributed k. The fluctuations in the population density 
is smaller on the BA model than on the ER graph, because the BA model has a broader 
degree distribution. 

To be more systematic, we compare population density fluctuation in the coexistence 
equilibrium. The fluctuation is measured by the standard deviation of the time series pi 
[see Fig. 11(b) and [1(c)] after transient, averaged over i = 0,1, and 2. Larger fluctuation 
means more unstable coexistence, and we examine how the size of the fluctuation depends 
on the amount of degree dispersion. In addition to the networks examined above, we use 
two types of networks that can create a range of degree dispersion. The first is the network 
with the two-point degree distribution. The standard deviation of the degree \J (k 2 ) — (k) 2 = 
\/p(l — p)\k\ — k 2 \. By varying ki/k 2 withp = 0.9 and (k) = 10 fixed, we can systematically 
create networks with a variety of y (k 2 ) — (k) 2 . The second is the network that has Gaussian 



Pk with mean (k), whose \J (k 2 ) — (k) can be also modulated. The results are summarized 
in Fig. H(d) for four types of networks (ER, BA, two-point, and Gaussian), excluding the 
R random graph because it does not sustain coexistence. Regardless of the network type, 
more dispersed degree distributions generally lead to more stable coexistence. 

In the mean field case, a smaller network with a stronger finite-size effect tends to drive 
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the population dynamics to the absorbing equilibrium where only one state survives 22]. 
The network effect on stability of coexistence is more manifested in this regime. In Fig. [21 
we show survival probabilities for some networks with n = 200, where the survival is defined 
by existence of all the three states. Figure [2] is consistent with Fig. QJd); coexistence is 
sustained for a longer period on networks with larger degree dispersion. 

III. MAY-LEONARD DYNAMICS ON NETWORKS 
A. Model and equilibrium 

Since neutrally stable oscillations of the RSP model may be sin gula r phenomena, we 



analyze another competition model proposed by May and Leonard [2l|]. The ML model 
represents dynamics of cyclically competing three species with natural death. Because of 
the natural death, vertices can take the vacant state. In a well-mixed population, the 
coexistence equilibrium and the periodic oscillation are both unstable. A trajectory of the 
population density approaches heteroclinic orbits on which at least one of the three species 
is extinct. Theoretically, one species is transiently and alternatively dominant with ever 
increasing periods in an infinite population. Practically, one species eventually wins due to 
the finite-size effect. 

As an interacting particle system, the ML model is a four-state system, with state 



representing the vacant site and 1, 2, and 3 representing cyclically dominating states 
The ML dynamics with heterogeneous contact rates are written as 



3- 
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(7) 



where po,fc = 1 — Pi,k ~ p2, k — P3,k- The first term in each equation indicates the rate at which 
a vacant site becomes colonized by state i (1 < i < 3). Supposing that a < 1 and (3 > 1, 
the second term and the third term of each equation represent the population increase and 
decrease due to the cyclic competition, respectively. For example, in the first equation, state 
1 outcompetes 2, whereas 1 is outcompeted by 3. Equation ([7]) has a coexistence equilibrium 
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given by 

P^ = e; = (a + /3 + l)" 1 (< = 1,2,3) (8) 

for all k. With p k = 5k,{k), the eigenvalues of the Jacobian matrix evaluated at (p* k , p 2k , 
P$u), disregarding a prefactor (a+fi+l)' 1 , are x = — (1+a+P), — l + (a+/3) /2±^/3(a — (3)i/2 
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. When a + (3 = 2, the ML model is essentially the same as the RSP model, and the 
coexistence is neutrally stable [x = —3, ±y/3(a — (3)i/2\. When a + (3 > 2 and a < 1 (or 
(3 < 1), the coexistence equilibrium is an unstable spiral, and the trajectory tends to a 
homoclinic orbit. 



B. Stability of coexistence equilibrium 

To investigate the network effect, we again consider the two-point degree distribution 
Pk = pSkM + (1 -p)fato- The Jacobian matrix at (p{ M , p{ M , p* M , p* M , p* 3ki , p% )k2 )\ where t 
denotes the transpose, is a block circulant matrix. Accordingly, the eigenmode has the form 
(v\, t>2, pvi, pv2, p 2 Vi, p 2 V2) t , v i, t>2 G C, where p is any solution to p 3 = 1. The corresponding 
characteristic equation is reduced to 

x 2 + ({(3 + pa + p 2 )(h + k 2 ) - (1 - p 2 ){f3 - x 



(k) 

+ (a 2 + f3 2 -af3-a-f3 + \)p 2 k x k 2 = 0. (9) 

For p — 1, Eq. ([9]) is a real equation, and two eigenvalues have the same real part that is 
equal to — (a + (3 + l)(ki + k 2 )/2 < 0. Because an eigenvalue of Eq. © for p = exp(27ri/3) is 
conjugate of one for p = exp(47n/3), it suffices to set p = exp(27ri/3). The larger real part of 
the solution to Eq. (jH)) is plotted in Fig. [3] for various p and < ki/k 2 < 1. The coexistence 
equilibrium is stabilized with negative real parts of the eigenvalues. This occurs when 
k\jk 2 = 0.1 and p = 0.9, 0.95. In this situation, a small number [= n(l — p)] of hubs have a 
large degree (= k 2 ) in comparison with most vertices with degree k\. This is reminiscent of 
long-tail pk typical of the scale-free networks. Excess heterogeneity {k\/k 2 < 0.05) destroys 
coexistence. In this situation, a majority of vertices with degree k\ are effectively isolated, 
and the network is close to the mean field case, or the R random graph with k = k 2 . 
When p < 0.5, the heterogeneity does not cause stability irrespective of k\jk 2 . This is 
because the contribution of the smaller subpopulation (proportion p) with the smaller degree 



S 



k = ki(< k 2 ) to dynamics is marginal, which again results in effectively homogeneous 
networks with k = k 2 . 



C. Numerical results 



For numerical simulations, we note that, in Eq. ([7]), a vacant site (state 0) is replaced 
by state % (1 < % < 3) at a rate of n^. A vertex in state 1 (2, 3) kills a neighboring state 2 
(3, 1) at a rate of f3 — 1. Then, the neighboring vertex is colonized by state 1 (2, 3) with 
probability (1 — a)/(/3 — 1) and turns empty (state 0) otherwise [121 ]. 

Dynamics for the R and ER random graphs with n = 5000 and (k) = 10 are shown 
in Figs. 11(a) andH^b), respectively. Because the stability condition for the ML model is 
more severe than that for the RSP model, one state shortly overwhelms the others on the 
ER as well as R random graph. However, the transients for the ER random graph, whose 
degrees are more dispersed than the R random graph, are longer. The BA model and 
the networks with two-point pt with parameters realizing the stable Jacobian matrix (but 
not the networks with Gaussian p/,) yield coexistence. Similar to the RSP dynamics, the 
amount of degree dispersion is strongly correlated with the amount of density fluctuation 
and stability of coexistence, irrespective of the network type [Fig. 11(c)] . 



IV. DISCUSSION 



A. Summary of the results 

We have examined population dynamics with cyclic dominance relationships on networks. 
The steady population density is independent of degree distributions of networks. However, 
stability of coexistence equilibria and dynamics depend considerably on networks. Het- 
erogeneity in degree distributions facilitates stable coexistence of different phenotypes. As 
touched upon in Sec. IH coexistence of competing species is desirable in, for example, eco- 
logical communities (biodiversity) and evolutionary games (survival of altruistic players). 
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B. Relations to synchronization of coupled oscillators 



The present mechanism of coexistence is related to synchronization of coupled oscillators. 
With spread degree distributions, each p^k evolves at a speed proportional to k, and p^i, 
Pi t 2, . . . are coupled by a sort of mean field feedback Oj. Then, the population dynamics are 
analogous to those of an ensemble of phase oscillators coupled by mean field interaction. In 
theory, coupled phase oscillators become desynchronized when the intrinsic frequency of the 
oscillators has a broad distribution relative to the coupling strength 24|. Oscillators with 
heterogeneous intrinsic frequencies correspond to vertices with heterogeneous k. In terms of 
competition dynamics on networks, asynchrony corresponds to stable coexistence of species 
where synchronous oscillations (large fluctuations in time) of the population density is sup- 
pressed. Desynchronization is known to suppress neutrally stable or unstable oscillations in 
ecological models with patchy populations, heterogeneous birth rates, and weak aggregation 
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Heterogeneity in degree distributions serves to stable coexistence via desynchronization 
even without other kinds of heterogeneity. The correspondence between asynchrony and 
coexistence may be exported to more general models, particularly to ones showing oscil- 
lations in well-mixed populations; oscillatory population densities are reminiscent of cyclic 
competition. 



C. Difference from spatial mechanisms of coexistence 



The scenario to coexistence unraveled here is distinct from those based on spatial struc- 
ture, heterogeneous environments, or heterogeneous individuals. In patchy habitats with 
heterogeneous environments or small diffusion 6], and in spatial structure with limited 

n n n □ 

[9|, 11Q|, llll, 112], multiple species can coexist by forming locally high densities of 



diffusion 



conspecifics in different subspaces [HJ, lUl, Il2j . This is the spatial mechanism of coexistence. 
In networks with dispersed degrees, multiple species can coexist on a network in a mixed 
manner without spatial segregation. 

Real networks of contacts are equipped with the clustering property, as is the case for the 
regular lattices and small-world networks, and high clustering elicits the spatial mechanism 
of coexistence. In addition, many networks own broad degree distributions represented by 
the scale-free distributions |l3l. |23||. The mechanism proposed in this work can cooperate 
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with the spatial mechanism to promote stability of coexistence. 



D. Difference from contagion dynamics 

There are many possible rules for interacting particle systems. In contagion processes, 
such as the percolation, the susceptible-infected-recovered model, and the contact process 
(susceptible-infected-susceptible model), degree dispersion affects dynamical aspects by, for 



example, accelerating disease propagation in initial stages 25| . More fundamentally, how- 
ever, epidemic thresholds (critical infection rates) are proportional to (k) / (k 2 ). Then, 
disease propagation on a global scale is more likely to occur on networks with more het- 
erogeneous contact rates with small (k) / (k 2 ) than on networks with rather homogeneous 
degrees such as the regular lattices and the random graph 4j, [23( . In contagion dynamics, 
the network influences the stationary state in addition to the dynamics, which contrasts to 
our results for the competition dynamics. 

Generally speaking, the positions of equilibria move when we cannot neglect at least one 
vpe of state-transition event whose occurrence rate is independent of neighbors' states n« 



20| | . Examples are spontaneous recovery and mutation. By contrast, the equilibria are 
invariant if all the state transitions are controlled by the neighbors' states, as is the case for 
the RSP model, the ML model, and also the voter model. However, the degree distribution 
does influence the stability of coexistence and hence the whole population dynamics. Our 
results are generalized to other population dynamics in which the rates of spontaneous 
transitions can be ignored. 
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Figure captions 

Figure 1: RSP dynamics on networks with n = 5000 and (k) = 10. The initial condition 
is given by the Bernoulli distribution with p — 0.7 and p\ = p 2 = 0.15, where pi is the 
proportion of vertices that take state i. The densities po (thin lines), p\ (moderate lines), and 
P2 (thick lines) are shown for (a) the R random graph. For (b) the ER random graph and (c) 
the BA model of the same size, only po is shown for clarity, (d) Fluctuation of population 
density as a function of the standard deviation of the vertex degree (ER, triangle; BA, 
horizontal line [({k 2 ) — (k) 2 ) 1 ^ 2 = 157.8]; Gaussian p k , crosses; two-point p k , circles). The 
variance of pi from time 150 through 300 averaged over i — 0, 1, and 2 defines the density 
fluctuation. 

Figure 2: Survival probabilities for the R random graph, the ER random graph, the BA 
model, and the networks with Gaussian pk with standard deviation 2 and 4 [corresponding 
to the crosses marked by arrows in Fig. [H(d)]. We set n = 200 and (k) = 10 for all the 
networks. The survival probabilities are calculated based on 1000 runs. 

Figure 3: Stability of the coexistence solution of the ML model. Real parts of the largest 
eigenvalues of the Jacobian matrix obtained from Eq. @ with a = 2/3 and (3 = 2 are 
presented for p = 0.1 (thinnest line), p = 0.3, p = 0.5, p = 0.7, p = 0.9, and p = 0.95 
(thickest line). We set (k) = 1 for normalization. 

Figure 4: ML dynamics on networks with n = 5000, (k) = 10, a = 2/3, and (3 = 2. The 
initial condition is given by po = 0, p\ = p<i = 0.25, and pz = 0.5. The R random graph 
(a) and the ER random graph (b) do not allow stable coexistence (po, dotted lines; pi, thin 
solid lines; p 2 , moderate solid lines; p 3 , thick solid lines), (c) Fluctuation of population 
density (BA, horizontal line [((k 2 ) — (k) 2 ) 1 ^ 2 = 157.8]; two-point pt, circles), defined by the 
variance of p; from time 150 through 300 averaged over i = 1,2, and 3. 
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